##### code chunk number 1: prepare###############################
options(digits=4)   #小数点位数
options(scipen=10)
graphics.off()      # 关闭图形设备
rm(list=ls())       # 清除对象
library(Cairo)      # R 图形输出pdf格式中文支持
library(mgcv)
library(splines)
library(lattice)
library(Matrix)
library(lme4)
library(nlme)
# END ###########################################################


##### 地区定义 ###################################################
huabei <- c("北京","天津","河北","山西","内蒙古")
dongbei <- c("辽宁","吉林","黑龙江")
huadong <- c("上海","江苏","浙江","安徽","福建","江西","山东")
zhongnan <- c("河南","湖北","湖南","广东","广西","海南")
xinan <- c("重庆","四川","贵州","云南","西藏")
#xinan <- xinan[-1]
xibei <- c("陕西","甘肃","青海","宁夏","新疆")
province <- c(huabei, dongbei, huadong, zhongnan, xinan, xibei)
whole <- c("全国",province)

expand.grid(1949:2011, whole)
# END ###########################################################

##### data import ###############################################
path <- "D:/Works/读博项目/Dissertation/Data/"
pdataquality <- read.csv(paste(path, "data1949.csv",sep=""))
pdataquality$gdp.rate <- pdataquality$gdp.rate/100 # 原数据没有包含%号，处理之
year.length <- 2011-1948
labor1989 <- pdataquality[pdataquality$year==1989,]$labor

#### 定义levels的顺序
pdataquality$district <- factor(pdataquality$district[drop=TRUE],
                                    levels=whole)

##### 1952 不变价GDP#

#海南1952-1978年GDP增长速度用广东省代替
pdataquality[pdataquality$district=="海南",]$gdp.rate[1:30] <-
    pdataquality[pdataquality$district=="广东",]$gdp.rate[1:30]

gdp.rate1952 <- unlist(
                tapply(pdataquality[pdataquality$year > 1952,]$gdp.rate,
       pdataquality[pdataquality$year > 1952,]$district,cumprod),
                use.names=FALSE)

gdp1952 <- ( rep(pdataquality[pdataquality$year==1952,]$no.gdp,
                 each=(2011-1952)) * gdp.rate1952)

no.gdp1952 <- pdataquality[pdataquality$year==1952,]$no.gdp

gdp1952 <- as.vector(rbind(matrix(nrow=3, ncol=32),no.gdp1952,
      matrix(gdp1952, nrow=2011-1952, byrow=F), deparse.level=0))

pdataquality <- cbind(pdataquality, gdp1952)

#### 海南缺失1978年之前数据处理
# 广东平减指数
gd.cum.rate1953.2011 <- cumprod(pdataquality[pdataquality$district=="广东",]$gdp.rate[-(1:4)])
gd.cum.rate1949.1952 <- rev(1/cumprod(pdataquality[pdataquality$district=="广东",]$gdp.rate[(4:2)]))
gd.cum.rate1952 <- c(gd.cum.rate1949.1952, 1, gd.cum.rate1953.2011)

gd.gdp1952 <- (pdataquality[pdataquality$district=="广东"
                           & pdataquality$year==1952,]$no.gdp*
               gd.cum.rate1952)

gd.deflator <- (pdataquality[pdataquality$district=="广东",]$no.gdp[-1]/
                (pdataquality[pdataquality$district=="广东",]$no.gdp[-63]*
                 pdataquality[pdataquality$district=="广东",]$gdp.rate[-1]))
gd.cum.deflator1953.2011 <- cumprod(gd.deflator[-(1:3)])
gd.cum.deflator1950.1951 <- rev(1/cumprod(gd.deflator[3:2]))
gd.cum.deflator1952 <- c(gd.cum.deflator1950.1951, 1, gd.cum.deflator1953.2011)

#海南1979-2011 环比 平减指数
hn.deflator1979.2011 <- (pdataquality[pdataquality$district=="海南" &
                             pdataquality$year>=1979,]$no.gdp/
                (pdataquality[pdataquality$district=="海南"&
                             pdataquality$year>=1978,]$no.gdp[-length(1978:2011)]*
                 pdataquality[pdataquality$district=="海南"&
                             pdataquality$year>=1979,]$gdp.rate))

#海南1950-1978 环比 平减指数 用广东省数据代替
hn.deflator <- c(gd.deflator[1:length(1950:1978)], hn.deflator1979.2011)

#海南1950-2011 定基 平减指数
hn.cum.deflator1953.2011 <- cumprod(hn.deflator[-(1:3)])
hn.cum.deflator1950.1951 <- rev(1/cumprod(hn.deflator[3:2]))
hn.cum.deflator1952 <- c(hn.cum.deflator1950.1951, 1, hn.cum.deflator1953.2011)

#海南1952年不变价GDP：只有1978-2011数据
hn.gdp1952 <- c(NA,pdataquality[pdataquality$district=="海南",]$no.gdp[-1] /
    hn.cum.deflator1952)

pdataquality[pdataquality$district=="海南",]$gdp1952 <-hn.gdp1952

# END ###########################################################

#### 1990 不变价资本存量K
k.deflator1990 <- pdataquality$k.deflator1952/
rep(pdataquality[pdataquality$year==1990,]$k.deflator1952, each=year.length)
no.stock.k <- pdataquality$stock.k1952*pdataquality$k.deflator1952
stock.k1990 <- no.stock.k/k.deflator1990
pdataquality <- cbind(pdataquality, stock.k1990)


pdataquality1990 <- pdataquality[pdataquality$year %in% 1990:2010,]#1990-2010
tapply(pdataquality$no.gdp, factor(pdataquality$year), sum, na.rm=T)


#### 1990 不变价GDP
gdprate1990 <- pdataquality1990[pdataquality1990$year==1990,]$gdp.rate
pdataquality1990[pdataquality1990$year==1990,]$gdp.rate <- 1
cumgdprate1990 <- stack(tapply(pdataquality1990$gdp.rate, pdataquality1990$district, cumprod))[,1]
gdprate1990 -> pdataquality1990[pdataquality1990$year==1990,]$gdp.rate
gdp1990 <- (rep(pdataquality1990[pdataquality1990$year==1990,]$no.gdp, each=21)*cumgdprate1990)
pdataquality1990 <- cbind(pdataquality1990, cumgdprate1990, gdp1990)

# END ###########################################################

##### 对重庆的处理，基本原则：合并到四川省#############################

#### 名义gdp直接相加
pdataquality1990[pdataquality1990$district=="四川", ]$no.gdp<-
    (pdataquality1990[pdataquality1990$district=="四川", ]$no.gdp
+pdataquality1990[pdataquality1990$district=="重庆",]$no.gdp)

#### 实际gdp直接相加
pdataquality1990[pdataquality1990$district=="四川", ]$gdp1990 <-
    (pdataquality1990[pdataquality1990$district=="四川", ]$gdp1990
+pdataquality1990[pdataquality1990$district=="重庆",]$gdp1990)

#### 受教育程度平均加权
pdataquality1990[pdataquality1990$district=="四川"
                 & pdataquality1990$year>1996, ]$edu <-
    (pdataquality1990[pdataquality1990$district=="四川"
                      & pdataquality1990$year>1996, ]$edu
     *pdataquality1990[pdataquality1990$district=="四川"
                      & pdataquality1990$year>1996, ]$labor
+pdataquality1990[pdataquality1990$district=="重庆"
                  & pdataquality1990$year>1996, ]$edu*
     pdataquality1990[pdataquality1990$district=="重庆"
                  & pdataquality1990$year>1996, ]$labor)/
 (pdataquality1990[pdataquality1990$district=="四川"
                      & pdataquality1990$year>1996, ]$labor
+pdataquality1990[pdataquality1990$district=="重庆"
                  & pdataquality1990$year>1996, ]$labor)

#### 劳动力直接相加
pdataquality1990[pdataquality1990$district=="四川"
                 & pdataquality1990$year>1996, ]$labor <-
    (pdataquality1990[pdataquality1990$district=="四川"
                      & pdataquality1990$year>1996, ]$labor
+pdataquality1990[pdataquality1990$district=="重庆"
                  & pdataquality1990$year>1996, ]$labor)

#去掉重庆
pdataquality1990 <- pdataquality1990[pdataquality1990$district != "重庆",]
pdataquality1990$district <- factor(pdataquality1990$district[drop=TRUE])

# END ###########################################################


##### 劳动力为年初和年底的平均 ######################################
labor.lst <- tapply(pdataquality1990$labor,
                   pdataquality1990$district, identity)
labor.df1990.2010<- data.frame(matrix(unlist(labor.lst), nrow=21))
labor.df1989.2009<- rbind(na.exclude(labor1989), labor.df1990.2010[-21,])
labor <- unlist((labor.df1989.2009 + labor.df1990.2010)/2)
pdataquality1990$labor <- labor

country <- pdataquality1990[1:21,]  # 全国独立出来
pdataquality1990 <- pdataquality1990[-(1:21),] # 各地区放在一起
# END ###########################################################


#### 去掉多余的levels
pdataquality1990$district <- factor(pdataquality1990$district[drop=TRUE])

pdq <- pdataquality1990
names(pdq)[9] <- "L"
names(pdq)[10] <- "EDU"
names(pdq)[12] <- "K"
names(pdq)[14] <- "GDP"
## DATA PREPARE END #############################################



##### Plot ######################################################

#### 1952-2011 名义GDP 地区总和与全国 差异
district.no.gdp.sum <-
    tapply(pdataquality[pdataquality$district != "全国",]$no.gdp,
           factor(pdataquality$year[-(1:63)]), sum, na.rm=T)
no.gdp.dif <- (district.no.gdp.sum
               - pdataquality[pdataquality$district=="全国",]$no.gdp)/
    (pdataquality[pdataquality$district=="全国",]$no.gdp)

#### 1952-2011 1952不变价GDP 地区总和与全国 差异
district.gdp1952.sum <-
    tapply(pdataquality[pdataquality$district != "全国",]$gdp1952,
           factor(pdataquality$year[-(1:63)]), sum, na.rm=T)
gdp1952.dif <- (district.gdp1952.sum
               - pdataquality[pdataquality$district=="全国",]$gdp1952)/
    (pdataquality[pdataquality$district=="全国",]$gdp1952)


## 差异作图

CairoPDF(file="chap6-fig-1.pdf", width=5.5, height=4.8)
plot(x=1952:2011,y=gdp1952.dif[-(1:3)], type="l", col="blue",
     xlab="年份", ylab="GDP", lty=4, family="AdobeSongStd-Light")
lines(x=1952:2011, y=no.gdp.dif[-(1:3)],col="red", lty=1)
text(x=c(1960,1960),y=c(0.3,0.265),c("名义GDP", "1952年不变价GDP"),
     adj=c(0,0),
     family="AdobeSongStd-Light")
legend(x=1950, y=0.35, legend=c("", ""),
       lty=c(1,4), col=c("red", "blue"), bty="n",
       )
# abline(h=0, lty=2, col="black")
# plot(pdq[c("GDP", "K", "L", "EDU")])
dev.off()

plot(x=1952:2011,y=gdp1952.dif[-(1:3)], type="l", col="blue")
lines(x=1952:2011, y=no.gdp.dif[-(1:3)],col="green", lty=4)
abline(h=0, lty=2, col="red")

plot(pdq[c("GDP", "K", "L", "EDU")], pch=".")

#### dotplot
set.seed(1234543)
print(dotplot(reorder(district, no.gdp) ~ no.gdp, pdq,
              ylab = "省份", jitter.y = T, pch = 21,
              xlab = "GDP(亿元)",
              type = c("p", "a"),
              family="AdobeSongStd-Light"
              ))

pdf(file="chap6-gdp-dotplot.pdf", width=7, height=5, family="GB1")
set.seed(1234543)
dotplot(reorder(district, no.gdp) ~ no.gdp, pdq,
              ylab = "省份", jitter.y = T, pch = 21,
              xlab = "GDP(亿元)",
              type = c("p", "a"),
#              family="AdobeSongStd-Light"
              )
dev.off()

set.seed(1234543)
print(dotplot(reorder(district, GDP) ~ GDP, data=pdq,
              ylab = "省份", jitter.y = F, pch = 21,
              xlab = "GDP(亿元)",
              type = c("p", "a")))
# END ###########################################################

##### Model #####################################################

#### Pooled Linear Model
lm.orig.1 <- lm(GDP ~ 1 + K + L + EDU, data=pdq)
summary(lm.orig.1)
AIC(lm.orig.1)

op <- par(mfrow = c(2,2), mar=c(4,4,2,2))
plot(x=fitted.values(lm.orig.1),y=resid(lm.orig.1),
     xlab="拟合值", ylab="普通残差")
plot(x=pdq$K, y=resid(lm.orig.1),
     xlab="K", ylab="普通残差")
plot(x=pdq$L, y=resid(lm.orig.1),
     xlab="L", ylab="普通残差")
plot(x=pdq$EDU, y=resid(lm.orig.1),
     xlab="EDU", ylab="普通残差")
par(op)

lm.orig.2 <- lm(GDP ~ 1 + K + L*EDU, data=pdq)
summary(lm.orig.2)
AIC(lm.orig.2)

lm.orig.3 <- lm(GDP ~ 1 + K + L:EDU, data=pdq)
summary(lm.orig.3)
AIC(lm.orig.3)

anova(lm.orig.1, lm.orig.2, lm.orig.3)

plot(lm.orig.2, which=c(1), add.smooth=FALSE)

op <- par(mfrow = c(2,2), mar=c(4,4,2,2))
plot(x=pdq$K, y=resid(lm.orig.2))
plot(x=pdq$L, y=resid(lm.orig.2))
plot(x=pdq$EDU, y=resid(lm.orig.2))
plot(x=(pdq$L)*(pdq$EDU), y=resid(lm.orig.2))
par(op)

#### Variance Structure
lm.orig.2 <- gls(GDP ~ 1 + K + L*EDU, data=pdq)

vf.Fixed <- varFixed(~K)

vf.Ident <- varIdent(form=~1|district)

vf.Power <- varComb(varPower(form=~K),
                    varFixed(~K)
#                    varPower(form=~L),
#                    varPower(form=~EDU)
                    )

vf.Exp <- varComb(varExp(form=~K),
                  varFixed(~K)
#                    varExp(form=~L),
#                    varExp(form=~EDU)
                  )

vf.ConstPower <- varComb(varConstPower(form=~K),
                    varConstPower(form=~L),
                    varConstPower(form=~EDU))

gls.orig <- gls(GDP ~ 1 + K + L+EDU, data=pdq,
                weights=vf.Ident) # 经过试验Power方差函数效果最好
summary(gls.orig)

plot(gls.orig,which=c(1), add.smooth=FALSE)


graphics.off()
CairoPDF(file="chap6-fig-gls-resid.pdf", width=5.5, height=5)
op <- par(mfrow = c(2,2), mar=c(4,4,2,2))
plot(x=fitted.values(gls.orig),y=resid(gls.orig, type="normalized"),
     xlab="拟合值", ylab="标准化残差",
     cex=0.5,
     family="AdobeSongStd-Light")
plot(x=pdq$K, y=resid(gls.orig, type="normalized"),
     xlab="K", ylab="标准化残差",
     cex=0.5,
     family="AdobeSongStd-Light")
plot(x=pdq$L, y=resid(gls.orig, type="normalized"),
     xlab="L", ylab="标准化残差",
     cex=0.5,
     family="AdobeSongStd-Light")
plot(x=pdq$EDU, y=resid(gls.orig, type="normalized"),
     xlab="EDU", ylab="标准化残差",
     cex=0.5,
     family="AdobeSongStd-Light")
par(op)

coplot(resid(gls.orig)~K|district, data=pdq)

#### Log Transformed Data
lm.trans.1 <- lm(log(GDP) ~ 1 + log(K) + log(L) + EDU,
                 data=pdq)
summary(lm.trans.1)
AIC(lm.trans.1)

lm.trans.2 <- lm(log(GDP) ~ 1 + log(K)
                 + log(L) + EDU + log(L*EDU), data=pdq)
summary(lm.trans.2)
AIC(lm.trans.2)

lm.trans.3 <- lm(log(GDP) ~ 1 + log(K) + log(L*EDU),
                 data=pdq)
summary(lm.trans.3)
AIC(lm.trans.3)

anova(lm.trans.1, lm.trans.2, lm.trans.3)

plot(lm.trans.2, which=c(1))


#### Mixed Model
## lme4 Mixed-effects modeling with R
(fm1 <- lmer(GDP ~ 1 + (1|district),
             data=pdq, REML=TRUE))

(fm2 <- lmer(GDP ~ 1 + (1|district) + (1|year),
             data=pdq, REML=TRUE))

dev.off()
qrr2 <- dotplot(ranef(fm2, postVar = TRUE), strip = FALSE)
print(qrr2[[1]], pos = c(0,0,1,0.6), more = TRUE)
print(qrr2[[2]], pos = c(0,0.4,1,1), more = TRUE)

dev.off()
print(xyplot(GDP ~ year | district, data=pdq,
             layout = c(6,5), type = c("g", "p", "r"),
             xlim = c(1990,2010), ylim = c(0, 15000),
             xlab = "",
             ylab = "GDP"))

## Random Intercept
fixed.1 <- gls(GDP~K+L+EDU+district, data=pdq, weights=vf.Ident)
summary(fixed.1)

gls.table <- xtable(summary(fixed.1)$tTable[1:4,],
                   caption="",
                   label="",
                   )
names(gls.table) <- c("估计值","标准误","$t$值", "$p$值")
print(gls.table,floating=F,
      caption.placement="top",
      table.environment=NULL,
      booktabs=T,
      sanitize.text.function = function(x){x}, # 数学公式
      only.contents=T, # 控制是否有table和表格环境
      )

mixed.intercept.1 <- lme(GDP~K+L+EDU,data=pdq,
                         random=~1|district,
                         weights=vf.Ident,
                         )
summary(mixed.intercept.1)


mixed <- lme(GDP~K+L+EDU, random=~0+K+EDU|district, data=pdq,
                       weights=vf.Ident,
                       )
summary(mixed)

gls.table <- xtable(summary(mixed)$tTable,
                   caption="",
                   label="",
                   )
names(gls.table) <- c("估计值","标准误","自由度", "$t$值", "$p$值")
print(gls.table,floating=F,
      caption.placement="top",
      table.environment=NULL,
      booktabs=T,
      sanitize.text.function = function(x){x}, # 数学公式
      only.contents=T, # 控制是否有table和表格环境
      )

varcov.table <- xtable(getVarCov(mixed)[1:2,1:2], digits=4)
K.d <- sqrt(getVarCov(mixed)[1,1])   # K的随机效应的标准差
EDU.d <- sqrt(getVarCov(mixed)[2,2]) # EDU的随机效应的标准差
K.EDU.cor <- getVarCov(mixed)[1,2]/(K.d*EDU.d)

align(varcov.table) <- c("","","")
print(varcov.table, floating=F,
      caption.placement="top",
      tabular.environment="bmatrix",
      booktabs=T,
      sanitize.text.function = function(x){x}, # 数学公式
      only.contents=F, # 控制是否有table和表格环境
      hline.after=NULL, include.rownames=FALSE, include.colnames=FALSE,
      )

mixed.intercept.2 <- lme(GDP~K+L*EDU,
                        random=~1|district, data=pdq)
summary(mixed.intercept.2)

mixed.intercept.3 <- lme(GDP~K+L:EDU,
                        random=~1|district, data=pdq)
summary(mixed.intercept.3)

anova(mixed.intercept.1, mixed.intercept.2, mixed.intercept.3)


## The Marginal Model等价于Random Intercept模型
gls.1 <- gls(GDP~K+L+EDU,
             correlation=corCompSymm(form=~1|district), data=pdq)
# K, L,EDU同时加入会导致收敛问题，故去掉L
summary(gls.1)

gls.2 <- gls(GDP~K+L*EDU,
             correlation=corCompSymm(form=~1|district), data=pdq)
summary(gls.2)

gls.3 <- gls(GDP~K+L:EDU,
             correlation=corCompSymm(form=~1|district), data=pdq)
# K, L:EDU同时加入会导致收敛问题，故去掉L:EDU
summary(gls.3)

anova(gls.1, gls.2, gls.3)


## Random Intercept and Random Slope
mixed.1 <- lme(GDP~K+L+EDU,
                        random=~1+K+EDU|district, data=pdq)
# K, L,EDU同时加入会导致收敛问题，故去掉L
summary(mixed.1)

mixed.2 <- lme(GDP~K+L*EDU,
                        random=~1+K+L*EDU|district, data=pdq)
#只有这个模型系数均为正
summary(mixed.2)

mixed.3 <- lme(GDP~K+L:EDU,
                        random=~1+K+L:EDU|district, data=pdq)
summary(mixed.3)

anova(mixed.1, mixed.2, mixed.3)


## The Marginal Model不等价于Random Intercept and Random Slope Model
gls.1 <- gls(GDP~K+L+EDU,
             correlation=corCompSymm(form=~1+K+EDU|district), data=pdq)
# K, L,EDU同时加入会导致收敛问题，故去掉L
summary(gls.1)

gls.2 <- gls(GDP~K+L*EDU,
             correlation=corCompSymm(form=~1+K+L*EDU|district), data=pdq)
summary(gls.2)

gls.3 <- gls(GDP~K+L:EDU,
             correlation=corCompSymm(form=~1+K|district), data=pdq)
# K, L:EDU同时加入会导致收敛问题，故去掉L:EDU
summary(gls.3)

anova(gls.1, gls.2, gls.3)
## END ###################################################


## Model Selection ############################################
lm.orig <- gls(GDP~K+L*EDU, method="REML", data=pdq)
mixed.intercept <- lme(GDP~K+L*EDU, method="REML",
                        random=~1|district, data=pdq)
mixed <- lme(GDP~K+L+EDU, method="REML",
                        random=~1+K+EDU|district, data=pdq)
AIC(lm.orig, mixed.intercept, mixed)
anova(lm.orig, mixed.intercept, mixed)
summary(mixed)

mixed.1 <- lme(GDP~K+L+EDU, method="REML",
                        random=~1|district, data=pdq)
summary(mixed.1)

mixed.2 <- lme(GDP~K+L+EDU, method="REML",
                        random=~1+K|district, data=pdq)
summary(mixed.2)

mixed.3 <- lme(GDP~K+L+EDU, method="REML",
                        random=~1+L|district, data=pdq)
summary(mixed.3)

mixed.4 <- lme(GDP~K+L+EDU, method="REML",
                        random=~1+EDU|district, data=pdq)
summary(mixed.4)

## mixed.5 不收敛
## mixed.5 <- lme(GDP~K+L+EDU, method="REML",
##                         random=~1+K+L|district, data=pdq)
## summary(mixed.5)

mixed.6 <- lme(GDP~K+L+EDU, method="REML",
               random=~1+K|district, data=pdq)
summary(mixed.6)

mixed.7 <- lme(GDP~K+L+EDU, method="REML",
                        random=~1+L+EDU|district, data=pdq)
summary(mixed.7)

## mixed.8 不收敛
## mixed.8 <- lme(GDP~K+L+EDU, method="REML",
##                         random=~1+K+L+EDU|district, data=pdq)
## summary(mixed.8)

anova(lm.orig, mixed.1, mixed.2, mixed.3, mixed.4, mixed.6, mixed.7)
## 最佳模型为 mixed.6
## END ###################################################


## 利用残差图检验是否存在异方差 ##############################
plot(mixed.6, which=c(1), add.smooth=FALSE)

op <- par(mfrow = c(2,2), mar=c(4,4,2,2))
plot(x=fitted.values(mixed.6), y=resid(mixed.6))
abline(h=0,lty=3)
plot(x=pdq$K, y=resid(mixed.6))
abline(h=0,lty=3)
plot(x=pdq$L, y=resid(mixed.6))
abline(h=0,lty=3)
plot(x=pdq$EDU, y=resid(mixed.6))
abline(h=0,lty=3)
par(op)

op <- par(mfrow = c(2,2), mar=c(4,4,2,2))
xyplot(resid(mixed.6, type="normalized")~fitted.values(mixed.6),
       panel=function(x,y){
    panel.grid(h=-1, v= 2)
    panel.points(x,y,col=1)
    panel.loess(x,y,span=0.5,col=1,lwd=2)})

xyplot(resid(mixed.6, type="normalized")~pdq$K,
       panel=function(x,y){
    panel.grid(h=-1, v= 2)
    panel.points(x,y,col=1)
    panel.loess(x,y,span=0.5,col=1,lwd=2)})

xyplot(resid(mixed.6, type="normalized")~pdq$L,
       panel=function(x,y){
    panel.grid(h=-1, v= 2)
    panel.points(x,y,col=1)
    panel.loess(x,y,span=0.5,col=1,lwd=2)})

xyplot(resid(mixed.6, type="normalized")~pdq$EDU,
       panel=function(x,y){
    panel.grid(h=-1, v= 2)
    panel.points(x,y,col=1)
    panel.loess(x,y,span=0.5,col=1,lwd=2)})
par(op)

## END ###################################################


## 加入方差结构 ############################################
vf.Ident <- varIdent(form=~1|district)

vf.Power <- varPower(form=~(L*EDU)) # K，L,EDU和积不收敛

vf.Exp <- varExp(form=~K)# K收敛，L,EDU和积不收敛,

vf.Fixed <- varFixed(~K)  # K收敛，L，EDU不收敛,积收敛

vf.ConstPower <- varConstPower(form=~K) # 三个变量均不收敛

mixed.var <- lme(GDP ~ 1 + K + L+EDU, method="REML",
                 random=~1+K+EDU|district,
                 data=pdq, weights=vf.Fixed) # 经过试验固定方差函数效果最好

summary(mixed.var)

plot(mixed.var,which=c(1), add.smooth=FALSE)


op <- par(mfrow = c(2,2), mar=c(4,4,2,2))
xyplot(resid(mixed.var, type="normalized")~fitted.values(mixed.var),
       panel=function(x,y){
    panel.grid(h=-1, v= 2)
    panel.points(x,y,col=1)
    panel.loess(x,y,span=0.5,col=1,lwd=2)})

xyplot(resid(mixed.var, type="normalized")~pdq$K,
       panel=function(x,y){
    panel.grid(h=-1, v= 2)
    panel.points(x,y,col=1)
    panel.loess(x,y,span=0.5,col=1,lwd=2)})

xyplot(resid(mixed.var, type="normalized")~pdq$L,
       panel=function(x,y){
    panel.grid(h=-1, v= 2)
    panel.points(x,y,col=1)
    panel.loess(x,y,span=0.5,col=1,lwd=2)})

xyplot(resid(mixed.var, type="normalized")~pdq$EDU,
       panel=function(x,y){
    panel.grid(h=-1, v= 2)
    panel.points(x,y,col=1)
    panel.loess(x,y,span=0.5,col=1,lwd=2)})
par(op)

## END ###################################################

#### GAMM
gamm.1 <- gamm(GDP ~ 1 + K + s(L)+EDU,
                 random=list(district=~1),#~1+K+EDU|district
                 data=pdq, weights=vf.Power) ##选择L为光滑样条为尝试结果

plot(gamm.1$gam)
plot(gamm.1$lme)
anova(gamm.1$gam)
anova(gamm.1$lme)
summary(gamm.1$gam)
summary(gamm.1$lme)
AIC(gamm.1$lme)

gamm.2 <- gamm(GDP ~ 1 + K + s(L)+EDU,
                 random=list(district=~K),#~1+K+EDU|district
                 data=pdq, weights=vf.Fixed)

plot(gamm.2$gam) # 此时L 为直线
plot(gamm.2$lme)
anova(gamm.2$gam)
anova(gamm.2$lme)
summary(gamm.2$gam)
summary(gamm.2$lme)
AIC(gamm.2$lme)


gamm.3 <- gamm(GDP ~ 1 + K + s(L)+EDU,
                 random=list(district=~1+K+EDU),#~1+K+EDU|district
                 data=pdq, weights=vf.Fixed) ## 不收敛, EDU也不收敛

plot(gamm.3$gam)
plot(gamm.3$lme)
anova(gamm.3$gam)
anova(gamm.3$lme)
summary(gamm.3$gam)
summary(gamm.3$lme)
AIC(gamm.3$lme)

#### Incoportate Temporal Correlation
xyplot(GDP~year|district,col=1,data=pdq)

gamm.4 <- gamm(GDP ~ 1 + K + s(L)+EDU,
               random=list(district=~1),
               data=pdq, weights=vf.Fixed,
               correlation=corAR1(form=~year|district))

gamm.5 <- gamm(GDP ~ 1 + K + s(L)+EDU,
               data=pdq, weights=vf.Fixed,
               random=list(district=~1),#~1+K+EDU|district
               correlation=corARMA(value=c(0.3,-0.1,0.2,0.4),
               form=~year|district, p=2, q=2)) ## 最佳p=2,q=2


gamm.6 <- gamm(GDP ~ 1 + K + s(L)+EDU,
               data=pdq, weights=vf.Fixed,
               correlation=corARMA(value=c(0.3,-0.1,0.2,0.4),
               form=~year|district, p=2, q=2)) ## 最佳p=2,q=2

gamm.7 <- gamm(GDP ~ 1 + K + s(L)+EDU,
               data=pdq, weights=vf.Fixed,
               random=list(district=~1+EDU),#~1+K+EDU|district
               correlation=corARMA(value=c(0.3,0.2,0.3,0.4),
               form=~year|district, p=2, q=2)) ## 最佳p=2,q=2


gamm.8 <- gamm(GDP ~ 1 + K + L+EDU,
               data=pdq, weights=vf.Fixed,
               random=list(district=~1+EDU),#~1+K+EDU|district
               correlation=corARMA(value=c(0.3,0.2,0.1,0.4),
               form=~year|district, p=2, q=2)) ## EDU变量不显著


mixed.var <- gamm(GDP ~ 1 + s(K) + L+EDU, data=pdq,
                  random=list(district=~1+K+EDU),
#                  weights=vf.Fixed,
#                  correlation=corARMA(value=c(0.4,0.2),
#                  form=~year|district, p=2, q=0)
                  ) # 经过试验固定方差函数效果最好


gamm.9<- gamm(GDP~K+s(L)+EDU, method="REML",
               random=list(district=~1+K),
               weights=vf.Fixed, data=pdq)
summary(gamm.9$lme)

gamm.9<- gamm(GDP~K+L+s(EDU), method="REML",
               random=list(district=~0+K),
               weights=vf.Ident, data=pdq)

summary(gamm.9$lme)
plot(gamm.9$gam)

AIC(gamm.1$lme, gamm.2$lme, gamm.4$lme, gamm.5$lme,
    gamm.6$lme, gamm.7$lme, gamm.8$lme, gamm.9$lme, gamm.10$lme, mixed.var)


plot(gamm.4$gam)
plot(gamm.4$lme)
anova(gamm.2$gam)
anova(gamm.2$lme)
summary(gamm.2$gam)
summary(gamm.4$lme)

plot(gamm.3$gam)


gamm.10 <- gamm(GDP~1+K+L+EDU, method="REML",data=pdq,
               random=list(district=~1+K),
#               weights=vf.Fixed,
#               correlation=corARMA(value=c(0.1,0.3),
#               form=~year|district, p=1, q=1)
                )

summary(gamm.10$lme)
plot(gamm.10$gam)
plot(gamm.10$lme)

#### 最佳模型
## 需要考虑自变量系数正负，是否显著，残差图形状
